//
// perc2 - calculation of the fractal dimension of correlated
// bond percolation cluster hulls
//
// Copyright (C) 2009, 2010 Indrek Mandre <indrek(at)mare.ee>
// http://www.mare.ee/indrek/perc2/, http://code.google.com/p/perc2/
// 
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
// 

#include <stdio.h>
#include <stdlib.h>

#include "perc2.hpp"

#define SQDIV 200
#define RDIV 720
#define SAMPLE 60

long double calc_cov (double x1, double y1, double x2, double y2)
{
  long double ret = 0;
  for ( int i = -1; i <= 1; i++ )
    {
      for ( int j = -1; j <= 1; j++ )
        {
          long double bcx = 0.5 + i;
          long double bcy = 0.5 + j;
          long double dist2 = (x1 - bcx) * (x1 - bcx) + (y1 - bcy) * (y1 - bcy);
          if ( dist2 > 1.5 * 1.5 )
              continue;
          long double v1 = v2dd(dist2hf (dd2v(dist2)));
          dist2 = (x2 - bcx) * (x2 - bcx) + (y2 - bcy) * (y2 - bcy);
          if ( dist2 > 1.5 * 1.5 )
              continue;
          long double v2 = v2dd(dist2hf (dd2v(dist2)));
          ret += v1 * v2;
        }
    }

  return ret;
}

int main ()
{
  FILE *fp = fopen ("calccellcov.csv", "w");
  for ( size_t li = 0; li <= SAMPLE; li++ )
    {
      //double len = 1.0;
      double len = li * (3.0 / SAMPLE);
      long double cum_cov = 0;
      for ( size_t yi = 0; yi < SQDIV; yi++ )
        {
          double y = 1.0 / ( 2 * SQDIV) + yi * (1.0 / SQDIV);
          for ( size_t xi = 0; xi < SQDIV; xi++ )
            {
              double x = 1.0 / ( 2 * SQDIV) + xi * (1.0 / SQDIV);
              for ( size_t ri = 0; ri < RDIV; ri++ )
                {
                  double phi = ri * (2 * M_PI / RDIV);
                  //printf ("%g %g %g\n", x, y, phi);

                  double x1 = x;
                  double y1 = y;
                  double x2 = x + len * cos(phi);
                  double y2 = y + len * sin(phi);
                  cum_cov += calc_cov (x1, y1, x2, y2);
                }
            }
        }
      printf ("%g,%.12Le\n", len, (RND_VARIANCE * cum_cov / (SQDIV * SQDIV * RDIV)));
      fprintf (fp, "%g,%.12Le\n", len, (RND_VARIANCE * cum_cov / (SQDIV * SQDIV * RDIV)));
    }
  fclose (fp);
  return 0;
}

